clear all; close all; clc
 addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/');
%----------- Load shocks using two identification schemes ----------------%
  NumberConsecutiveshocks = 2; % 1 previous shocks that is positive, and the 2nd one is positive too 
FromM2Q = 0;  
cd ./DataShocks/
  LoadShocksMandQpositiveEPUrevisedMAC; sample=[ 1961 1; 2019  12];
cd ../  
   noEPUhistorical = 1;
if noEPUhistorical
                                       %sample=[ 1985 1; 2019  12];
                                        sample=[ 1986 1; 2019  12];
 %ss_ = find(dates_EPU ==198501);
  ss_ = find(dates_EPU ==198601);
  se_ = find(dates_EPU ==201912);
  uEPU   =   uEPU(ss_:se_);
  logEPU = logEPU(ss_:se_);
  IndPreviousNShocksEPU = IndPreviousNShocksEPU(ss_:se_);
  clear ss_ se_
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%--------------------------------------%
clear startsample endsample 

T = size(data,2);
 P  =  6; %                including in the LP regression P lags of all variables in the system
 
  plottype = 'riskfr';      ii=1; 
                           %ii=2; % for Wu-Xia shadow rate
% plottype = 'Output';      ii=3; %
% plottype = 'Inflat';     %ii=8; % 
%                           ii=9; % 
%                          %ii=4; % CPI
%                          %ii=5; % PCE
% plottype = 'SP500' ;     %ii=6; %      
%                           ii=7; % real SP500 
%% IRF 
h1 = 0 ;  %start LP at h=0 or h=1 (h=1 if impose no contemporaneous impact)
H  = 12*3; %#horizon of IR

% We are interested in the estimation of the dynamic multiplier of y_t+h with respect to a change in x_t for h ranging from 0 to H
%------------ IRF of GDP growth to a shock variable ----------------------%
 y  = data(:,ii);    %endogenous variable
% if ii == 8 || ii == 9
%     y  = [NaN(12,1); data(13:end,ii)-data(1:end-12,ii)];    %endogenous variable
% end
%-------------------------------------------------------------------------%
  Remove_Large_Shocks = 1;
    x  = uEPU;      
%/////////////////////////////////////////////////////////////////////////%
  w  = [lagmatrix( data(:,unique([ii 1 10    ])  ), 1:P )                       ];
 if ii == 9
  w  = [lagmatrix( data(:,unique([ii 6 1 10 3])  ), 1:P                        )];
 end
  w( ~isfinite(w) ) = 0;
OKforReg = isfinite(x) & isfinite(y);
for cc=1:size(w,2)
    OKforReg = OKforReg & isfinite(w(:,cc)) ;
end  
  x = x(OKforReg,:); 
  y = y(OKforReg,:); 
  w = w(OKforReg,:);
  
% This shrinks the estimated dynamic multiplier to a polynomial of order r - 1
r = 3; %(r-1)=order of the limit polynomial (so r=3 implies the IR is shrunk towards a quadratic polynomial)

lambda =  20; %value of penalty
%-------------------------------------------------------------------------%
% k-fold cross validation 
  k_fold=5;
%-------------------------------------------------------------------------%       
lambda =  20; %value of penalty 
 lp = locproj(y,x,w,h1,H,'reg');             %IR from (standard) Local Projection
slp = locproj(y,x,w,h1,H,'smooth',r,lambda); %IR from Smooth Local Projection

figure(41)
%%
 scalingPLOTS = 1.0;
%scalingPLOTS = 3.0; %Remember that BBD compute model-implied responses of industrial production and employment to a 90-point upward EPU innovation
  plot( 0:H,scalingPLOTS*[ lp.IR slp.IR],'LineWidth',3); hold on;  
  plot( 0:H,scalingPLOTS* slp.Interval_up,'r--','LineWidth',2); plot( 0:H ,scalingPLOTS*slp.Interval_down,'r--','LineWidth',2);
plot( 0:H , zeros(H+1,1) , '-k' , 'LineWidth' , 2 )
grid
xlim([0 H])

switch plottype
    case 'Output' 
        title('Industrial production','FontSize',16);     %set(gca,'XTick',t(4:4:end),'FontSize',12);
        ylim([-1.6 1.0]);set(gca,'YTick',[-1.6:0.5:1.6],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'riskfr'
        title('Risk-free rate','FontSize',16);%set(gca,'XTick',t(4:4:end),'FontSize',12);
        %ylim([-0.4 0.4]);set(gca,'YTick',[-0.6:0.2:0.6],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'SP500'
        title('Stock market','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
        ylim([-4.0 1.4]);set(gca,'YTick',[-4.0:0.5:1.4],'FontSize',12);ylabel('Percent','FontSize',12)
    case 'Inflat'
        title('Inflation','FontSize',16); %set(gca,'XTick',t(4:4:end),'FontSize',12);
        %ylim([-5.4 5.4]);set(gca,'YTick',[-5.2:0.5:5.2],'FontSize',12);ylabel('Percent','FontSize',12)    
    otherwise
        warning('Unexpected plot type. No plot created.')
end
legend('IR_{lp}','IR_{slp}',                             'Location','Best')

keep slp ii 
fname = sprintf('EPUBBD_LP_linear_ForReplic.mat');
if     ii==1
    slpShortRate= slp.IR; slpShortRateUP= slp.Interval_up; slpShortRateDown= slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown')    
elseif ii==3
    load(fname)
    slpOutput     = slp.IR; slpOutputUP     = slp.Interval_up; slpOutputDown     = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown')
elseif ii==6 || ii==7
    load(fname)
    slpStockMkt = slp.IR; slpStockMktUP = slp.Interval_up; slpStockMktDown = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown','slpStockMkt','slpStockMktUP','slpStockMktDown')
elseif ii==9
    load(fname)
    slpInflation = slp.IR; slpInflationUP = slp.Interval_up; slpInflationDown = slp.Interval_down; 
    save(fname,'slpShortRate','slpShortRateUP','slpShortRateDown','slpOutput','slpOutputUP','slpOutputDown','slpStockMkt','slpStockMktUP','slpStockMktDown','slpInflation','slpInflationUP','slpInflationDown');
end